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^ ; ABSTRACT 

■ Mature neutron stars are cold enough to contain a number of superfluid and super- 

conducting components. These systems are distinguished by the presence of additional 
O ■ dynamical degrees of freedom associated with superfluidity. In order to consider models 

' with mixtures of condensates we need to develop a multifluid description that accounts 

for the presence of rotational neutron vortices and magnetic proton fluxtubes. We also 
need to model the forces that impede the motion of vortices and fluxtubes, and under- 
stand how these forces act on the condensates. This paper concerns the development 
of such a model for the outer core of a neutron star, where superfluid neutrons co-exist 
with a type II proton superconductor and an electron gas. We discuss the hydrody- 
' namics of this system, focusing on the role of the entrainment effect, the magnetic 

fjS . field, the vortex/fluxtube tension and the dissipative mutual friction forces. Our final 

results can be directly applied to a number of interesting astrophysical scenarios, e.g. 
associated with neutron star oscillations or the evolution of the large scale magnetic 
field. 

O ' 



jj ■ 1 INTRODUCTION 

Neutron stars host the strongest magnetic fields known in nature. The magnetic field of a typical pulsar has an (already) 
enormous intensity of ~ 10 12 G, but this is nothing compared to magnetars which can reach, or even exceed, the extraordinary 



intensity of ~ 10 G (IDuncan fc Thompsonlll992D . Given this, it is not surprising that the magnetic field plays a key role 



in neutron star physics. As far as neutron star dynamics is concerned, the presence of the magnetic field is likely to affect 
large-scale phenomena like pulsar glitches and free precession. It may also alter the properties of the star's various oscillation 
modes. In fact, if the field is sufficiently strong, it may provide the main restoring force acting on the fluid elements, leading 
to Alfven-type oscillation m odes. Such effects m ay already be observed in quasi-periodic oscillations associated with the tails 



of magnetar giant flares (see IStrohmaverl l|2007l ) for a review) 



In its vast majority, existing w ork on the dynamics of mag netised neutron stars makes use of standard magnetohydro- 



dynamics (MHD) (see, for example. IShapiro fc Teukolskvl |l983|)) which is based on a "single-fluid" hydrodynamical model. 



Although this is a reasonable initial approach to the neutron star problem, it is of questionable validity in more realistic 
scenarios. This is easy to see if we recall that below a critical temperature threshold T c ~ 5 x 10 9 K (promptly reached after 
a neutron star's birth) the liquid matter in neutron star interiors is expected to undergo a phase transition to a quantum 
"ordered" state where the bulk of the neutrons (protons) become superfluid (superconducting) through the formation of 
Cooper pairs. As a result of this phase transition the macroscopic hydrodynamical behaviour of neutron star matter is that 
of a multifluid system with more than one distinct fluid degree of freedom. In addition, the presence of a superconductor in 
neutron star matter means that the properties of the macroscopic magnetic field are likely to differ from the familiar plasma 
physics behaviour. 

The number of distinct "fluids" appearing below T c depends on the composition of matter. Since the outer core of a 
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canonical neutron star is primarily composed of neutron-proton-electron (npe) matter, one can identify three fluidi^]. The 
key point is that superfluidity allows for relative flows among interpenetrating matter components. In the inner neutron star 
core the situation could be further complicated by the presence of exotic particle species like hyperons or deconfined quarks. 
However, even though this leads to very interesting problems, we will not consider the deep core in this paper. 

A consistent formalism for neutron star magnetohydrodynamics must be based on a multifluid model which accounts for 
superfluidity and superconductivity. The main aim of the present paper is to formulate such a m odel. Naturally, this subject 



has been addressed in the past, in particular in a series of papers by Mendell and Lindblom (|Mendell fc Lindbloml 11991 
Mendel] 1991 . 1998). Our analysis builds on their pioneering work. 



We consider a neutron star as a three fluid system (superfluid neutrons, superconducting protons and normal electrons) 
at ze ro temperature, coupled thro ugh entrainment, mutual friction forces and the magnetic field. Despite this common ground 
with iMendell fc Lindbloml i|1991T ) we present a substa ntial revision of the subject. In particular, we correct a conceptual 
error in the original calculation by Mendel] ( 1991 , 19981 ) , regarding the energy contribution of the so -called London field (the 
charac teristic magnetic field induced in a rotating superconductor). This error was first pointed out by Carter. Prix fc Langloisl 
(2000), although they did not attempt to reformulate the MHD problem. As we will show in this paper, this error has no 
serious repercussions for most of the final neutron star MHD equations. The reason for this is the intrinsic weakness (several 
orders of m agnitude smaller than the stellar magnetic field) of the London field. Hence, most of our final equations are similar 
to those of Mendel] ( 1991 . 19981 ). A notable exception is Ampere's law, which will be shown to depart significantly from its 
usual form as a result of the London field. 

The main new contribution of our work is the derivation of the magnetic force in the MHD equations for super- 
fluid/superconducting neutron stars. In particular, we provide a general proof that the standard electromagnetic Lorentz 
force is eliminated fro m these e quations, be ing replaced by a tension force due to the proton vortices (fluxtubes) . This result 
was also obtained by Mendel] ( 1991 . 199a), although it was demonstrated only for a plane- wave model. In addition to the 
derivation of the superconducting magnetic force, we provide a detailed discussion on the various vortex-mediated mutual 
friction forces that appear in the equations of motion. For each of these forces we provide an explicit (al beit phenomeno logical) 
expression and discuss its properties and expected st rength. This analysis builds on previous work by iMendelll (|l9911 ). 

The final difference between our work and that of Mendell ( 1991 . 1998h concerns the choice of variables. We d o not use the 
"orth odox" approach to superfluid hydrodynamics, where the superfluid "velocities" in fact represen t momenta dKhalatniko " 
1965|). Ins tead, we prefer to use the canonical formulation developed by Carter and collaborators ( Carter fc Langloisl Il99 



Prixll2004l : TAndersson fc Comejbood ) which makes a clear distinction between kinematical veloc ities and m omenta, avoiding 



conceptual confusion. Since the two approaches are mathematically equivalent (see, for example. iPrbd ((200J)), opting for one 
or the other is a matter of taste. Yet, we believe that the v ariational formulation offers a deeper insight into the physics of 
multifluid systems. In particular, it uses conserved fluxes, see Andersson fc Corner! d200S ) . This approach is also advantageous 
since it is readily extended to fully general relativistic models (jAndersson fc Comerll2007l ). 



2 AN OVERVIEW OF NEUTRON STAR SUPERCONDUCTIVITY 

To set the stage for the discussion of the MHD of superfluid and superconducting neutron star s, we first provide a qualitative 



over view of neutron star superconductivity. Our discussion draws heavily on the classic work of Bavm. Pethick fc Pines ( 1969h 



and lMendelll |l99lh . which in many ways still represents the state-of-the-art in this area. However, we also touch upon some 



modern aspects of superconductivity. 

The effect of superconductivity on neutron star dynamics is two-fold. Firstly, the protons decouple from the rest of 
the matter and provide a distinct macroscopic fluid degree of freedom, in a manner analogous to the superfluid neutrons. 
Secondly, the stellar magnetic field exhibits a "non-classical" behaviour due to the Meissner effect, the expulsion of magnetic 
flux from superconducting matter. In order to formulate the MHD equations for superconducting neutron stars it is necessary 
to understand the implications of these effects. 

One may expect to find many similarities between the well established physics of terrestrial superconductors and neutron 
stars. In reality, this is only partially the case. Laboratory systems exhibit a complete flux expulsion due to the Meissner effect 
for any magnetic field strength below a critical threshold. This is not the case for neutron stars, where the superconducting 
matter can be penetrated by any finite magnetic field. Basically, the magnetic field is metastable with the global Meissner 
effect o perating on a secular timescal e (of the order of Myrs), due to the high ambient electric conductivity of neutron star 
matter ( Bavm. Pethick fc Pinesl 1969f ) . 



The magnetic flux penetrates the matter by redistributing itself into filamentary structures within which proton super- 
conductivity is suppressed. The topology and properties of the resulting regions, occupied by normal protons, depend on the 



1 This is strictly true only if we ignore finite temperatu r e effects. Otherwise the entropy associated with finite temperature excitations 
may provide an additional degree of freedom (|Prixll2004j ; lAndersson fc Come r 2008, 2012). 
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type of superconductivity jTillev fe Tillevlll99d ; iTinkhamlfljij ). The familiar dichotomy between type I and type II super- 
conductivity also applies for neutron stars (at least when there is a single superconducting particle species). Which type of 
superconductivity prevails depends on the relative size of two characteristic lengthscales: the proton coherence length £ p , on 
which Cooper-pairing can be broken by proton quantum fluctuations, and the penetration length A* , on which the magnetic 
field is exponentially shielded inside superconducting matter. The type of superconductivity is determined by the ratio 



\/2A» 



(1) 



In the case of type II superconductivity re s < 1 and it is energetically favourable for the magnetic field to penetrate the 
superconductor by forming an array of quantised proton vortices[f|, each carrying a single flux quantum (j>o = hc/2e ~ 
2 x 10~ 7 Gem 2 . Locally, each proton vortex has a cylindrical geom etry with a nor mal proton core of radius £ p . The detailed 
structure of proton vortices in neutron stars has been discussed bv lMendelll |l99ll ). who finds 

1/3 /z P \ 1/3 A0 9 K 



£ p « 5 x ICT 



(2) 



and (see Appendix) 
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where p±4 = p/10 14 gcm -3 is the normalised mass density, x p is the proton fraction, m p is the bare proton mass and m* 
is the effective mass acquired by each proton as a result of the entrainment. Finally, T cp is the (density dependent) critical 

temperature for proton superconductivity. 

As in a laboratory superconductor, the upper critical field H C 2 is given by 1 Tinkham 199ot ). 
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This corresponds to a state where the proton vortices are densely packed to the point where their cores "touch" , leading to a 
global destruction of superconductivity. That is, we have 
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where d p is the proton vortex separation. 

In the laboratory, the vortex state is only accessible for fields in the range H c i < B < H C 2, where the lower critical field 
is given by 



Hd = 



4ttA 2 



. 'A. 
Tog I — 

sp 



Taking (1/2) log (A./| p ) w 1, cf. iTinkhaml l| 19961 ). we have 



2ttA 2 



: 4 X 10 1 



(6) 



(7) 



Each proton vortex is associated with a magnetic field Bp that points along the (local) direction of the vortex core and is 
distributed within a radius ~ A» . That is, we have 



9o 



H cl (8) 

A lab oratory superconductor can only be penetrated by an imposed magnetic field between the critical values ( Tillev fc Tillevl 
1990). As we have already mentioned, the neutron star problem is different. In this case, the vortex state is accessible also for 



B < H c i, although it is metastable in this region. Hence, the transition at H c \ is not likely to be as important for astrophysical 
systems. In addition, it is straightforward to show that type II superconductivity should be relevant for neutron stars. By 
combining the expressions for the coherence and penetration lengths we find, 
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For typical par ameters pertaining to the out er core of a mature neutron star we have k s < 1, i.e. type II superconductivity 
should prevail ( Bavm. Pethick fc Pines 1969). However, at the critical density 



Peril ~ 3 x 10 
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(10) 



2 Rather than using the terminology proton fiuxtubes (fiuxoids) and neutron vortices, we will refer to both as "vortices" associated with 
each fluid. This makes sense since, in the neutron star case, the neutron vortices are magnetised. Moreover, both sets follow from the 
quantization of the (canonical) momentum circulation. 
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Figure 1. This figure provides a superconductivity "cross section" of a mature neutron star core (assuming npe composition) as a 
function of the density. We plot the critical temperature T cp for proton condensation (normalised to 10 10 K), the critical fields H c 2, H c i 
(normalised to 10 16 G) for type II superconductivity, and the ratio k s which determines the type of superconductivity. Note that the 
vertical axis is dimensionless. The left (right) panel corresponds to a model of strong (weak) proton superfluidity (see text for details). 
The vertical dashed line indicates the p cr it density for the type II — >■ I transition. 



beyond which k s exceeds unity, a type II — > I phase transition should take place. This transition ought to be relevant for 
neutron stars, since the critical density can easily be reached in the inner part of the core, see Figure Q] 

In a type I region the formation of proton vortices is energetically prohibited. Instead, the magnetic field is distributed 
into regions of normal pr otons with l arge s patial dimensions (compared to the fluxtube cross-section) and without a universal 
shape (see for example, Sedrakian 1 20051 )). The actual shape of these regions is mainly determined by the magnetic field 
distribution when the protons condense (at the critical temperature T cp ). The lack of well-defined structures (like vortices) in 
a type I superconducting region is the main obstacle preventing the formulation of macroscopic MHD for this region. For this 
reason, we will not consider type I superconductivity further in this pape r. How e ver, s ince a portion of a neutron star core 
may contain type I matter, future work needs to address this problem (see I Jones (2006) for an initial discussion). 

So far, our survey of neutron star superconductivity has been somewhat idealised, and our development of the MHD 
equations will be at a similar level. This is natural, since this is our first investigation of the problem. However, it is useful to 
mention some other aspects of superconductivity which could be relevant (and perhaps of key importance) for more realistic 
neutron star models. 

One such aspect is the variation with density of the key physical parameters, like the length scales £ p and A*, the critical 
temperature T cp and the ratio n a . Using the relation A D ~ 1.75fcB?cp betwee n Tcp and the proton pairing energy gap A p and 
the approximate gaps discussed bv lAndersson. Comer fe Glampedakis l|2005l ) we can produce neutron star "superconductivity 
cross-sections" like that shown in Fig. [T] This figure spans a typical range between "strong" and "wea k" proton pairing (the 



chosen models correspond to the pairing cases "e" and "f" in I Andersson. Comer fe Glampeda kis (2005), right and left panels 



respectively), and illustrates the layered profile of proton superconductivity in a neutron star core composed of npe matter 
only. 

The results in Fig. [1] show that an attempt to model realistic neutron stars would have to account for the coexistence 
of type I, type II and normal regions and the detailed physics at the interfaces separating these regions. In addition, the 
figure shows that the critical field H C 2, as given by eqn. (UJ), can vary by up to an order of magnitude. While it may reach a 
maximum ~ 10 16 G it can drop to ~ 10 15 G as the transition density is approached. 

New features may also appear due to the coexistence of several mutually interacting superfluid and superconducting 
particle species. This possibility has usually been ignored in discussions of neutron star physics, for example, in the calculation 
of the structure of a vortex. Yet, it may turn out that these interacti ons are importan t . For t he "simple" two-component mixture 
of superconducting protons and superfluid neutrons recent work by Alford fe Good! (2008) suggests that the type II — > I phase 
transition is significantly affected by the presence of the neutron superfluid. The coupling between the two condensates results 
in a shift in the value of n s at which the transition takes place. Moreover, the coupling leads to the presence of a region 
populated by "higher order" proton vortices, each carrying a larger number of flux quanta. This result hints at the complexity 
of the real type II > I transition and the presence of an interface that could have impact on neutron star dynamics. 

The departure from standard superconductivity could be even more pronounced in the inner core, especially if it is 
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populated by exotic par ticles. The most li k ely sc enario allows for the presence of hyperons in superfluid/superconducting states 
(for a recent review see Schaffner-Bielicbl 1I2 OO8)). If we were to consider only the S~ and A hyperons, which are expected to 
appear first l Haensel. Potekhin fe Yakovlevl 2007 ). then we end up with a five- fluid system, containing two superfluids (n, A) 
two superconductors (p, S~) and the normal electrons. This is clearly a much more complicated multifluid system than that 
for the two-fluid outer core. Moreover, hyperon cores may exhibit the phenomenology of so-called "two-gap" superconductors. 
An example of such a system is liquid metallic hydro gen which exhibits a variety of vo rtex states, not observed in standard 
single-component superconductors (see, for example. iBabaev. Sudbo fc AshcrofiJ ((2004J)). This means that the high density 
part of the superconductivity cross-sections shown in Fig. [T]could be drastically different in neutron stars with hyperon cores. 

Another issue that has not yet been discussed in detail is th e character of superconduc tivity in magnetars. For magnetars, 
the exterior dipole magnetic field is estimated to be ~ 10 15 G ( Woods fe Thompsonll2004r ). The strength of the interior field 
is unkno wn but it should be comp arable, or even higher, than the exterior one, probably with a mixed poloidal/toroidal 
topology |Braithwaite!l2008 fe 200m . This is interesting because magnetic fields in the range ~ 10 15 — 10 16 G may be strong 
enough to suppress proton superconductivity (above H C 2)- Hence, it is not unreasonable to suggest that large portions of a 
magnetar core may not be superconducting. Another likely (and interesting!) possibility is that the interior magnetic field lies 
between the critical fields H c i and H C 2- Then the proton vortices will be so densely packed, in the sense that their spacing 
dp < A,, that they will interact. This regime of superconductivity has not received any attention at all in the context of 
neutron stars, but it is known fr om laboratory systems that m utual vortex interactions could lead to a new type of ultra-low 
frequency magnetic oscillations i De Gennes fc Matricon 19641 ). 



3 MULTIFLUID MHD WITHOUT VORTICES 



As a first step towards a formalism for mixed superfluids/superconductors, we discuss the derivation of the MHD equations 
for a self-gravitating multifluid system composed of a number of charged and neutral constituents (labeled by the indices 
{x, y}) coupled to an electromagnetic field. At this point, any number of constituents can be in a superfluid/superconducting 
state. Initially, we will not account for the presence of vortices. This will be done later, in Section [4] 

The discussion in this Section can be kept at an abstract level without the need to identify the various constituents. 
Of course, at the end of the day, the resulting MHD equations will be specialised for modelling the outer core of a neutron 
star composed of superfluid ne utron s, supe rconducting protons and normal electrons. It is worth noting that the m odel we 
will d iscuss , qhich is based on [ Prixl ( 2004 ) . represents the Newtonian limit of the relativistic model developed by Carter 



1 19891 ). see lAndersson fc Corned (|2007l ~ :or a review. By including charged constituents, it is also the natural extension of 



the flux-conservative model discussed by Andersson fc Comer 1 20061 ). However, we will not attempt to develop an analogous, 
generic, model for the dissipative problem. In fact, it is well-known that the development of such a model, accounting for the 
electromagnetic field, is a real challenge. 

The basic variables of the multifluid system are the number densities n x , the kinematic velocities v x , the mass m x and 
the charge g x carried by each particle of the x constituent. Given these, we can define mass and charge densities and the 
associated fluxes/currents according to, 



px = m x n x 
f x = g x n x , 



n x = n x i; x 
fx = ct x v x 



Summing over all components we get the total charge current, mass density and charge density, 

■ i \ ^ -i \ \ ^ 



(11) 
(12) 

(13) 



The system can be described by a Lagrangian C, which is a functional £(n x , n x , Aq, A 1 , $), where the gravitational potential 
is $, and Ao and A 1 represent the electromagnetic degrees of freedom (the vector potential). The latter parameters are related 
to the electric field Ei and the magnetic field Bi in the usual way, i.e., we have 



E, = ViA - -d t Ai, 
c 



/>- = > -' VjA, 



As a result, the two Maxwell equations 



are trivial identities. 



ViB' = 0, e^VjEk 



-d t B l 



(14) 



(15) 



Following IPrbd (|2005l ). we assume a Lagrangian with minimal coupling. This means that the "hydrodynamical" contri- 
bution can be separated from the electromagnetic one. In other words, we take as our Lagrangian 



C = £ H (n x , O + £v M {Ei, Bi) + £ grav + ( aA + ^-jiA x ) - p$ 



(16) 
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The electromagnetic and gravitational contributions to the Lagrangian are given by the usual expressions 

£em = ^(E 2 - B 2 ), £ grav = -^(V$) 2 (17) 
Meanwhile, the hydrodynamical piece is given by (|Prixll2004h 



£ H = ^ m ; n ' nx -g (is) 



2n x 

where the energy £ represents the "equation of state". In general this energy has a functional form £ — £(n x ,w 2 y ), where 
w xy = v * — v y is the relative velocity between the x and y constituents. 

The remaining two Maxwell equations are obtained from the variation of C with respect to Aq and A 1 . This leads to 

V,D l = 47R7, e l]k V 3 H k = —f + -d t D* (19) 

c c 

where the "dynamical" fields (D l ,H l ) are related to the "kinematical" fields (E l ,B l ) according to 

*r|| (20) 
Similarly, variation of C with respect to the gravitational potential $ leads to the Poisson equation 

V 2 $ = 4ttGp (21) 
The hydro dynamical equations of motion are derived by varying £ with respect to the basic fluid variables n x and n x (see 



Prix! ((2004J) for details). If the total particle number for each constituent is assumed fixed then one arrives at the particle and 



charge continuity equations, 

d t n x + V s n x = 0, d t a x + V*j£ = (22) 
For an isolated system, the conservation of momentum is then expressed by a collective Euler-type equation 



El e\ i r-ri X ilk X r-jl m I 

< n x dt7r x — n x V 7r — e rije&imV 7r x > 



(23) 



This equation features the canonical momenta 

dC x o x 



,;=^=rf + ^ (24) 



where we have introduced the hydrodynamical momenta (Prix 2004), 



<9n x 



^ x + ^e xy wf x j 



Pi = l^T = ^(". X + )' Exywf I (25) 



This expression demonstrates the well known fact that the entrainment effect between the different particle species, entering 
the formalism through the parameters 

' (26) 



xy p x \dw* 

can misalign the momentum p x with respect to the corresponding flux n x . 
The scalar function 7r x — dC/dn^ in the equations of motion is equal to 

1 2 



7r x = — n x — -m x v K + q*Ao — m x $ (27) 



where 

,,. = ( 

, dn 

are the usual chemical potentials. 

The momentum equation (|23[) can be re-written on an equivalent "force-balance" form (note that the signs here are 
convention) 



Hx ~ ^EMx + /3 X V 1 $ } = (29) 



If there are additional (external) forces acting on the system they will ap pear on the right-hand-side of this equation. The 
"hydrodynamical" and electromagnetic force densities in (|29|l are given by (|Prixll2005h , 



Fix = nxj (d t + « x V 3 )p x + V> x + mxJ2 (exywf) VX j 



(30) 
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fmx = ^ E % + ^"rii:, J (31) 

At this point it is worth noting that in this model the number of Euler equation-type terms in (|23[) or (|29[) is the same as 
the number of constituents. In this sense, each constituent can be identified with a "fluid". This identification can always 
be made in principle, but in many situations it is not appropriate. If two or more constituents are comoving (i.e. have same 
Hx) then it is obviously only their "conglomerate" which represents an identifiable fluid degree of freedom. This is the case 
when the dynamical timescale is much longer than the relaxation due to interparticle scattering or electromagnetic coupling. 
Superfluidity suppresses dissipative scattering and allows the components to have distinct dynamics. It is also worth pointing 
out that a single "constituent" fluid may have multiple dynamical degrees of freedom. This is the case for, for example, 4 He 
at low temperatures, where a superfluid condensate coexists with a gas of thermal excitations (phonons and rotons). The 
description of this system also requires two- fluid hydrodynamics (see Andersson fc Comet (2008) for a recent discussion). 

The model we have discussed is general in the sense that it applies to any multifluid sys tem interacting with an elec- 
tromagnetic field. When compared to the standard equations of plasma physics ( JacksorJ 19981 ). the formalism outlined here 
shows certain similarities but also significant differences. The main difference is d ue to the presence of entr ainment-induced 
couplings. This is known to be a generic property of interacting Fermi liquids (jAndreev fc Bashkinl 119751 ). Since one can 
associate the entrainment coefficient with an "effectiv e" mass for ea ch constituent, it is also a familiar effect from other 
mixing-problems (like gas bubbles moving in liquids, see Geurst ( 1988T )). Additional differences will emerge when constituents 
become superfluid/superconducting and the system is endowed with rotational and magnetic vortices. 



3.1 The MHD approximation 



The MHD approximation follows the assumption of macroscopic charge neutrality. This obviously requires the various charged 
con stituents to be c arriers of opp osite charge, and there must be an overall balance. As discussed in many standard textbooks 
(cf. Jackson ( 19981 ); Chen (1984)), charge neutrality holds to high precision at lengthscales larger than the so-called Debye 



length and for time intervals larger than the plasma period (~ 10 -21 s). Typically, these requirements are easily met when 
charged particles (like electrons) are sufficiently mobile that they neutralise any local charge. This renders the system oblivious 
to its true plasma origin. In addition to charge neutrality, the MHD approximation uses the simplifying assumption of "slow" 
fluid motion (compared to the speed of light) . It is easy to verify t hat, for the conditions i n neutron star interiors, the MHD 
approximation is justified at any scale relevant to hydrodynamics ( Easson fc Pethick 19791 ). 

Hereafter we will focus on a particular multifluid system: a three-constituent mixture of neutrons, protons and electrons, 
with respective labels x = {n, p,e}. This is the simplest model for the outer core of a neutron star. Generalisations to other 
situations would be, at least in principle, strai ghtforward. Consider, for example, the inclusion of muons. Then, due to the 
short timescale for electron- muon scattering (|Mendelllll99ll ). the electrons and muons will rapidly relax to be co-moving on 
average. It is thus easy to extend our mod el to cover this case, but we have opted not to do this (in order to keep the discussion 
clear). We refer the interested reader to iMendelll (|199ll ) for details and briefly remark, where appropriate, how our results 
would change if muons were present. 

Taking the charges to be q p — — q a — e and q n = 0, and imposing overall neutrality a = <r c + <r p = we must have 
Vi.f = and, 

n c = n p =>• f — eripWp C (32) 

This leads to Ampere's 



(33) 



The slow motion approximation then allows the elimination of the displacement current from (|19p 
law, 

Let us now consider the momentum equations for this system. Without any loss of generality, equation (|29|l can be separated 
by prescribing a coupling force F^ f between the neutrons and the combined charged components. At this point this force is 
unspecified, but we will later identify it with the vortex-mediated mutual friction force. This leads to the system 



(dt + v J n \7j)p' n + Vfin + mne n M)j n V'v° + m n V l $ > = F^t 



(34) 



where 



n p \ (d t + v J p\7 3 )pl + VV P - mpEpW^V 1 ^ \ + p c (d t + v 3 c Vj) v\ + n e VV e + (p P + Pe)V ! $ = -F^ t + F£ (35) 



(36) 



Fii 1 ilm ■ n 

L = -e ]iB n 
c 



is the standard Lorentz force density. 

In equations (|34[) - (135l) we have only considered the entrainment between the neutrons and the protons. The coefficients 
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£n = £n P and e p = £ pn are related since p n £ n = £ P p P , c.f. (|26|) . We assume a vanishing entrainment for the electrons, i.e., 
let £ex = 0, in order to simplify the problem. There is no compelling physical argument not to make this assumption. The 
situation is different for the entrainment between the baryons, which is due to the strong interaction. 

Even though the system contains three distinct constituents, it is possible to "reduce" it to an effective two-fluid system 
by exploiting the smallness of the electron inertia (compared to that of protons and neutrons). Mathematically, this amounts 
to neglecting the electron inertia in (|35(l . Thus, we obtain 



where we have introduced the combined proton-electron chemical potential, 

fi = Mp+Mc (38) 

Ttlp 

In the following we will use m = m n = m p . We arrived at (|37p by neglecting terms which are a factor m c /m <C I 
smaller than those we retain. In other words, these equations are accurat^f] up to about one part in 10 3 . The resulting Euler 
equation (|37[) is effectively a single-fluid equation, without explicit dependence on the electron velocity. The electron velocity 
enters only the total current j* which, however, can be eliminated by means of (|33[) . If we ignore the Lorentz force, e.g. 
if the protons and electrons are considered to be comoving, then the pair of equations (|34p and (|37[) reduce to the non- 
magnetic two-fluid model that has been used in a number of studies of superfluid neutr o n stars dynamics (see for example 
|PrixNl999l):IPrix fc Rieutordl ll 19991): I Andersson fc Comerl d200lll:lchamel fc Carterl JpOOrj^lGlampedakis. Andersson fc Jonesl 
l|2008 fc 2009h : lGlampedakis fc Anderssonl (|2009F 7 iPassamonti. Haskell fc Anderssonl |2009l )). 



3.2 The induction equation 



In order to complete the set of MHD equations, we need to provide a relation between the electric field and the fluid variables. 
Given such a relation, we can eliminate the electric field from Faraday's law (|15p and derive a final equation relating the 
magnetic field to the fluid variables, the so-called induction equation. 

We take the Euler equation for the electron fluid as our starting point. This equation can be separated from the collective 
equation of motion ([29} by introducing (at this point unspecified) coupling forces Fp B and F„ e between the electrons and the 
proton and neutron fluids, respectively. This means that we have 



Udt+viVj)vi + V' (Ac + <*>)] +■ 



c 



vjBk 



F' 



■ F 1 

- 1 pe 



Again ignoring the electron inertia we obtain an expression for the electric field 



E 1 



1 

-— € 

C 



1 3 ' 



v C jB k 



-V* (/i e + m c $) — ( 



Fl 



i 1 pc 



Inserting this result in (|15ft and eliminating the current, we get 



d t B l = e 



B r ' 



+ 



Pp 



where 



e 

mc 



The electron velocity can be eliminated via 



(39) 



(40) 



(41) 



(42) 



(43) 



en p p 47ra p p p 

Equation (|41|l is the ge neral expression for the induc tion equation in our multifluid systenQ. It differs from the standard MHD 
induction equation, c.f. Shapiro fc Teukolskv 1 19831 ). in two ways. Most apparently, the total coupling force F* = Fp C + F^ c 
does not necessarily take the form of a frictional force F c ! oc f (which in standard MHD would lead to an Ohmic dissipation 
term in (|4ip ). This point will become obvious when we discuss the nature of Fl in Section [5] In addition, we will later 
establish that H r (entering in (|41[) implicitly through the use of (|43[) ~) differs considerably from the standard result. In fact, in 
Section |4]3] we will establish that H l = b^, where &l is the London field. Both these effects may lead to important differences 
regarding, for example, the magnetic field evolution. 



3 In the presence of muons eqn. d 371) is still valid provided we redefine eqn. (38) according to 

PpVi/i ft! (pp + p c + p M )Vi/i = PpVi/ip + pcVi/ic + p^Vip M 

The equations are in this case accurate up to corrections of the order of (p e + Pn)/Pp ~ 10 -2 if we assume that n M < 0.1n p which is 
reasonable in the outer core of a neutron star. 

4 In fact, eqn. l|41|l is valid even if muons are present, but the forces F* e and F^ e should of course be reinterpreted as being due to 
interactions with the combined electron-muon fluid. 
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4 MHD IN THE PRESENCE OF VORTICES 

In order to model the conditions expected in the outer core of a mature neutron star, we need a formalism that accounts 
for neutron superfluidity and proton superconductivity. This is of key importance since fundamental features like the bulk 
rotation and the magnetic flux may be "locked" into a population of quantised neutron and proton vortices. These need to 
be incorporated in the formalism. 

The presence of Fermi-scale sized vortices leads to a physical system with multi-scale dynamics. However, as in the case 
of superfluid Helium, it is usually sufficient to consider the smooth-averaged properties of the neutron and proton vortex 
arrays. This is at least the case as long as one is mainly interested in large-scale dynamics. Averaging results in macroscopic 
quantisation conditions, the so-called Onsager-Feynman relations (equations (|44l) and (|45[) below), which involve the fluid 
velocities and the magnetic field. In the case of neutron stars it is obvious that one can only make progress by working at 
the averaged (macroscopic) levejf). There is no way that one will ever be able to study the dynamics of the immense number 
of individual vortices. Neither would one want to, for the result would likely be too complex to interpret. The averaging 
is essential, but comes at a price. The averaged model is essentially phenomenological, and the link to the key microscopic 
parameters is not always clear. In fact, one of the main challenges concerns the inference of macroscopic parameters from 
detailed microscopic models. 

In addition to the quantisation conditions, the presence of vortices leads to new forces in the hydrodynamical equations. 
Loosely speaking, these forces tend to be of two different kinds: (i) resistive "drag" forces resulting from the interaction 
between the vortices and the fluids, leading to effective "mutual friction" coupling terms in the Euler equations, (ii) forces 
due to the intrinsic "tension" of the vortices themselves. So far, we have ignored all such forces. The following Sections are 
devoted to the discussion of the forces and the associated modification of the MHD equations from Section [3] 

Before proceeding to the detailed analysis, it is important to make a comment on notation. In the rest of the paper, the 
various parameters entering the MHD formalism refer to averaged quantities and as such the adopted notation should formally 
distinguish them from the corresponding parameters of the preceding vortex-free formalism. However, in order to keep the 
presentation as simple as possible, we opt to use the same notation. In doing so, there should be no danger of confusion since 
the remainder of the paper is focussed on vortex-averaged MHD. 



4.1 The macroscopic quantisation conditions 

Due to the fundamental quantum nature of the superfluid/superconducting condensates, neutron and proton vortices will be 
quantised in such a way that they carry a single quantum of momentum circulation^. We will assume that both neutron and 
proton vortices are locally arranged in rectilinear arrayfl directed along unit vectors and k' p (we use a 'hat' to denote 
unit vectors). The two arrays have su rface densities M n and J\f p , respectively. At the hydrodynamics level, we then have 
( Mendell fc Lindblom|[l99ll ; |Prixll2005l ) , 

W i = ! £ «* v jP l=M a j n (44) 
m 

and 

w * = i e «»V j5 rJ = ^£ tf *V,-pS + = Af P n P (45) 

We have used k]. = kkJ, with k = h/2m. It is important to note that the quantised "vorticities" refer to the circulation of 
the (averaged) canonical momenta 7r* = p], + (q x /c)A l rather than the circulation of velocity. It is the canonical momentum 
which is related to the gradient of each condensate's wavefunction phase </> x , leading to the Onsager-Feynman quantisation 
condition 

jnidh = (h/2) j(V l ^)dh = ft/2 . (46) 

The vortices interact with their environment mainly through their intrinsic magnetic fields. These fields are of the order 
of 10 15 G and extend only a very short distance away from the vortex core (see Appendix). The magnetic field of the 
neutron vortices is generated by the current of protons that are entrained by the local circulation associated with the vortex 



This is in contrast with laboratory systems like cold atom condensates which often require an approach that "resolves" individual 
vortices. 

6 It would be relatively straightforward to extend our analysis to let each proton vortex carry a larger flux. However, we expect that the 
most common case will be that of single quantum vortices. This case should be energetically favoured in most situations. 

7 Although it is usually assumed that the vortices in a neutron star core are "aligned" , there is no reason why this should be so. In fact, 
from the many studies of turbulence in superfluid laboratory systems one may expect that the generic situation involves vortex tangles. 
The implications of such t angles for neutron star dynamics is not yet well understood. Ne vertheless, it is clear from the first efforts aimed 
at exploring this problem (Pcralta ct al. 2005 & 2006; Andcrsson, Sidcrv & Comer 2007) that turbulence could play a key role in many 
situations. 
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( Alpar, Langer fe Saulslll984h . The protons can react to flows on the very short lengthscale associated with the vortex core 
since they themselves form a condensate with a short coherence length. The mean-free path of the electrons is much larger, 
which means that they cannot respond collectively on this short lengthscale. Inferring the mean-free path (A e ) from results 
for electron-electron scattering shear viscosity (with coefficient rj), we have 

A c « -2- (47) 
since the electrons are expected to be relativistic. Scaling to typical parameter values, we then have 

A^xlO-^;' Vl^W^cm (48) 



,10 15 g/cm 3 7 V T 

This result implies that the electrons do not act as a "fluid" on the mesoscopic scale of the individual vortices. As a result, 
the entrainment between neutrons and protons leads to the neutron vortices being magnetized. 

The total averaged magnetic field B l in (|45p is a sum of three distinct parts: the averaged fields contributed by the 
neutron and pro ton vortex arrays and the so-called London field &£• The London field is a fundamental property of a rotating 
superconductor ( Tillev fc Tillev 1990T I. In contrast to a neutral superfluid which mimics bulk rotation by forming a vortex 
array, a superconductor can rotate without forming vortices. As we will soon see, the proton fluxtubes themselves are, in fact, 
not related to macroscopic rotation at all. This crucial difference between superfluids and superconductors is encoded in the 
Onsager-Feynman equations ((44)) and l|45[l. 

The total averaged magnetic field can be written 

B l = bl + B; + Bl (49) 
The averaged vortex fields are expressed in terms of the densities A/" x and the magnetic flux associated with a single vortex, 

Bl = A/"x<Mi (50) 
Each proton vortex carries a single quantum of magnetic flux 

he 

<Pp = 4>a = — (51) 

Meanwhile, the flux associated with a neutron vortex can vary continuously, as it is induced by the proton current that is 
entra ined by the neutron circulation. As discussed in the Appendix, each neutron vortex carries a flux ( Alpar. Langer fc Saulsl 



19841 ) 

4hx = - T ^ s —4>o (52) 

1 En 

where the minus sign appears because the magnetic field of a neutron vortex is antialigned with K n . Although the fluxes <p p 
and <f> n are comparable, the total magnetic field is entirely dominated by the proton vortex contribution since Af P exceeds Mb. 
by several orders of magnitude for typical neutron star parameters. 

By combining (|49[) with (|44|l and (|45|l we find that the London field is given by 

6 l = -J- (e ijk V jP l - -IP-e^Vjpl) (53) 



If the entrainment is assumed uniform, this simplifies to 

* ijk 



En 



where we have defined, for later convenience, 



K = ~^i< (54) 

(2r> 



(55) 



It is useful t o note that the entrainment par ameter can be expressed in terms of the effective proton mass m p . Since s 
1 — m*/m p ( Prix. Comer fc Andersson 20021 ) we find that 

— p p + prn*/m p ^ m* 
p — 2p p + p p m*/ m p m p 
This should be a good approximation in a neutron star core, where the proton fraction is small. 
In the case of uniform proton rotation with angular frequency Q p , eqn. (J54J) becomes 



(56) 



bl = -^ifit (57) 



p 



This result clearly reveals the link between the London field and the rotation of the proton superconductor. Inserting typical 
numerical values in (|57p shows that the London field is exceedingly weak compared to the typical magnetic field strengths in 
neutron stars, 



-i 



bL „ _p „ o.Ol ( — ] G (58) 
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4.2 A variational derivation 

The MHD equations from Section 13.11 need to be modified in the presence of neutron and proton vortices. The new set of 
equations must be cognizant of the quantisation conditions (|44[) and (|45[) and of the forces associated with the vortices on a 
mesoscopic scale. At the level of a single vortex a Magnus force is present when there is relative motion with respect to the 
macroscopic fluid flow. A tension force appears when the vortex "bends" , much like the tension in a piano wire. There is also 
a Lorentz force due to the interaction of the intrinsic vortex magnetic field with the charged fluids. 

Once the above pieces are combined in a consistent way, the end product is a set of macroscopic Euler equations featuring 
averaged force densities t x due to the vortices. That is, we have 

f>Jy (dt + «£Vj) K + + V 1 {fin + #) + e n «4 n VX J = F*a + t (59) 

and 

PvUdt + upVj) K " £ P «V) + V 1 (/i + *) - eptojnVwJ J = + Fl+t l p (60) 

The inclusion of these forces is, however, not the only necessary modification to the equations of motion. As we discuss later, 
cf. Section 14.41 there are also significant changes to Ampere's law and the Lorentz force F£ (due to the averaged macroscopic 
current j % interacting with the averaged total magnetic field B l ) for a type II superconductor. 

The "complete" equations of motion can be derived by a suitable extension of the L agrangian framework described in 



Section [3] This ap proach is closely related to t h e original variat i onal m ultifluid formalism of Be karevich fc Khalatnikovl (|196lh 



section |3| this approach is closely related to the original variational multinuic 
(later adopted bv lMendell fc Lindbloml Jl99lh: iMendelJ jl99li . Il99sh ) which 
and to the more recent framework of ICarter fe Langloisl (|1995l ). 



is based on the use of the total energy density, 



The presence of neutron and proton vortices in the multifluid system can be taken into account by adding an extra piece 
C v to the total Lagrangian, 

C = £ H + CEM + £grav + C v + (<tA + \jiA*\ - (61) 

This vortex term is a functional C v — £ v (n x , W£) of the canonical momentum vorticities and the number densities. Its specific 
form is given by the total averaged vortex contribution to the energy, which is obtained by simply taking the product of the 
energy per unit length £ vx of a single vortex and the corresponding vortex density 

£v = -^A/"x£vx (62) 

X 

The addition of C v is not the only required modification to the Lagrangian. When we discussed the variational approach to 
the magnetised fluid problem in Section[3]we made a point of separating the hydrodynamical and electromagnetic contributions 
to the Lagrangian. When vortices are present this separation is no longer appropriate. Since both neutron and proton vortices 
are magnetised, the piece (|62[) also contains the averaged energy due to the local magnetic fields. We need to modify £bm 
accordingly and account for the contribution of the "free" electromagnetic field, i.e. the London field. Thus, we have 

£em = {E 2 - bl) (63) 

Since, as we will shortly see, the London field is itself a functional ?>L(n x , B l , W x ) it is more natural to consider the combination 

b 2 

Co = C V -^- (64) 

07T 



when taking the variation of C. Then the total Lagrangian i^f] 



E 2 ( 1 \ 

L = £ H + £ g rav + — + crAo + -jiA 1 - p$ + Co = £ vf + C a (65) 
Sir \ c J 

where C v f represents all the pieces that retain the same form as in the vortex- free system, see Section [3] In particular, it is 
important to point out that the minimum coupling term in (|65[1 refers to the total averaged electromagnetic field, as given by 
eqn. 1(49)) . 

Before proceeding with the variational calculation it should be noted that, in (|61[) and (|65[) we have not accounted 
for contributions due to interactions between neighbouring proton and neutron vortices. Neutron vortex-vortex interac- 
tions originate from the long-range radial profile of the neutron circulation. The macroscopic manifestation is an elastic 
response of the vortex array to perturbations that do not involve vortex "bending". This effective elasticity leads to the 

8 In terms of the energy density used in the formalism of iMendell fc Lindbloml il99lT ); iMendelll il99ll , | 1998) the present modifications 
to the total Lagrangian amount to an energy expression So = b^/Six + A/x£vx, omitting the contributions that do not play a role in 
the present discussion. 



12 Glampedakis, Andersson & Samuelsson 



so-called Tkachenko waves. These ultra-low fre quency waves have recently been observed in laboratory superfluid systems 
i Coddington et allbooj iMizushima et al.l 12004 ). Whether they are relevant in the context of neutron stars is still an open 



issue, see lNoronha fe Sedrakianl ([2008) for a recent discussion. The main interaction between proton vortices, which is due to 
their short-range magnetic fields, could become important if the vortices "touch", in the sense that d p < A*. This is possible 
provided the magnetic field is sufficiently strong to satisfy H c i < B < H C 2- This regime could be relevant for magnetars with 
superconducting cores, but we will not consider it further here. 

Returning to the Lagrangian C and working out the variation we get 



5C = 5£. vi - {a?<W x + Cx<5n x | + -^-HiSB 1 



with the definitions 



As before 



A? 



d£ 
dWi 



dC 



Hi 



. dC A 8C 

47T-r-^- = 47T- 



(66) 



(67) 



(68) 



dB i dB i 
since the only independent piece of £ is Co- 

From (|66[) we can see that the chemical potentials will have additional electromagnetic contributions such that 

— > + Cx (69) 

In order to arrive at the specific form of the vortex forc es in the eq uations of motion we need to express the variations SW^ 
in terms of the Lagrangian fluid displacement field £ x IjPrbd 12004! ) . The required relation follows from the conservation of 
vorticity, which at the local level requires 



AW* = -W: 



6Wl = e ijk e klm V j (&W?) 



(70) 



where A represents a Lagrangian variation. Given this result, it is straightforward to show that the vorticity terms manifest 
themselves as the forces 



l W^V i A^ = W^(V i Af-V I A x ) 



in the equations of motion. This agrees with the conclusions of Mendelll (1991). 

It is relatively easy to see why the force should take this form. We can rewrite the original fluid equation (|30[) as 

ijk x 



/Hx = "xftPx + TlxV 1 ( jU x - ^"x + «xPj 



K, X x-rl m 

Vj EfclmV p x 



(71) 



(72) 



Recalling that the variational analysis |Prixj[2004j) used a Lagrangian Cn{nl.), where n x = n x t)J are the individual fluxes, we 
rewrite the last term as 

ijk xt— tI m /f7(i\ 

— e tkimrijV p x (73) 
Comparing this to (I71[) we see that the latter has exactly the form that one might expect from the variational analysis. 



With the vortex terms £ x and t\ accounted for, the canonical equations of motion (|23[) become (after combining eqns. (145 
459j>, © and CD}), 



V*« - Cx) - e tJk {n 3 + e 3ps V p \:} e fcIm V l 7r™ } = 



(74) 



Note that the scalar function -kq is now defined as ttq = dC v i/dn^ and, therefore, it is still given by (127[) , 

In order to proceed to the calculation of the explicit form of Co and the resulting vortex terms t x an d Cx we need to 
ensure that has the required functional form. This is easily achieved by combining the conditions (144[) and (I45[) , 



B l 



K ^ 



Using this result, together with (recall that J\f x — W x /k) 



Co = --Vw„-ij : 

K z — ' Sty 



(75) 



(76) 



9 This is analogous to the loca l mass conservation law for each species, which takes the form ((Friedman & Schut3 Il97 
lAndersson, Comer fc Grosar3l2004l 'l 

An x = <5ra x + £ x Vjn x = -n x (Vjf x ) . 
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we find 



K = — w x 

K 



9o 



where in the last equation we have defined 



47ra p £ vx 



4tt£ v 



K 4>o 

We have also also used VV X = k x . Then from eqns. (|71[) and (|67|) it follows that 

1 



t 



47ra D 



-WZ 



v J (H vx wr)-v 1 ( J ff vx w J x ) 



00 



+ Vi 



¥>x ,L 



and 



Cx = 



4- 9px Wy + dp, 

y^x 



1 - En 



(77) 



(78) 



(79) 



(80) 



The parameters H vx play the role of effective magnetic fields. In particular, from the definition and the explicit results for 
the vortex energies (eqns. Ij82p and (|83[) below) we can identify H vp — H c i and estimate H vn ~2x 10 15 pi4 G. These results 
suggest that H vn /H c \ ~ 10. 

O ur fin al results (|79[) and (|80p for the vortex forces t], and potentials £ x do not agree with those obtained by iMendell 



1 199ll . 1 19981 ). The difference is, however, limited to te rms which depend on the gene rally weak London field. As we have 
already pointed out, based on the original discussion by Carter. Prix fc Langlois ( 2000h ■ the origin of this discrepancy can be 
traced back to an erroneous vortex energy contribution used by Mendell. Translated into our notation this energy (£ vx )m is 
written as 



(£vx)m = £vx + -r-Kib 

47T 



X / 2 

L 



(81) 



From a conceptual point of view, the appearance of a macroscopic variable like in the expression f or the mesoscopic vortex 
energy is suspicious. As we explicitly show in Appendix [A] (generalising the original calculation of Carter. Prix fc Langloisl 
(|2000h ) it is, indeed, incorrect. The correct expressions for the neutron and proton vortex energies are given by eqns. (|A28[) - 
(fA30T> . 

For an isolated neutron vortex the energy is well approximated by (see eqn. (|A30|) ) 



'I 



(82) 



47T 1 - £„ 

Recall that the required energy integral is cut off both at the inner edge (letting the hollow vortex core be represented by the 
coherence length £ n ) and at a distance where the assumption that the vortex is isolated breaks down (taking this to be the 
vortex separation <i n ). Similarly, the energy of a proton vortex is 

'A,' 



47T£* 



lot! 



(83) 



where the energy integral has a cut-off at a distance comparable to the London penetration length A* 



4.3 Ampere's law for a type II superconductor 

In normal matter (or in vacuum) the dynamical field H l is typically identified with the magnetic field itself, i.e. H l = B l . 
However, this is no longer true for a type II superconductor. From the definition (|68p we obtain the, possibly surprising, result 

^ = -47rg => H* = l& (84) 



As far as we are aware, ICarter. Prix fc Langloisl ifcood ) were the first to arrive at this result, pointing out that the usual 



identification H l = B r should not be expected to hold in type II superconductors. In retrospect, this result makes perfect 
sense given that the bulk of the magnetic field is carried by the quantized proton fluxtubes. The system is analogous to a 
magnetized medium rather than a magnetic classical fluid. 

It follows that Ampere's law for a type II superconductor is now given by (neglecting friction forces) 

e VA = — j = -w P e (85) 

c c 

This is obviously different from the standard resullF°l. 

The result (|85[l allows us to draw some general conclusions. First, the presence of the weak London field b]^ instead of the 

10 Our form for Ampere's law does not agree with the corresponding result of Mendell (1991, 1998) who makes the usual identification 
H l = B 1 , as a result of using the erroneous vortex energy 1811 1. 
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full field B' in (|85[1 suggests that the total macroscopic (averaged) current will be significantly weaker than what we would 
expect in normal matter, given that in neutron stars we would expect 6l >C B. 

Second, in the case where the proton fluid is uniformly rotating, and the entrainment parameter e* is also assumed 
uniform, the combination of (|85[) and (|57[) leads to the somewhat surprising conclusion that the proton and electron fluids are 
comoving, v p — v z c . In the more realistic scenario where e* is allowed to vary inside the star (while preserving axisymmetry 
with respect to the spin axis) we obtain 

i4 e «2A^fip^V (86) 

where r is the cylindrical radial coordinate. 



4.4 The total magnetic forces in the equations of motion 

We now return to the complete equations (|59l) and (|60[) and write them as 



and 



- -I- 

mf i r m; 



i p 



having defined £ x = Cx/ m and the total effective magnetic force 

F^ ag = t p + Fi- P P V% 

exerted on the charged fluids. 

Using the definition of the Lorentz force 

F i = l e ^ jlBm = _L B i ( V H i _ v 

C 4-7T V 

we make the identification l|84|) and obtain 

F i = ^b 3 (v.foL-v 1 ; 

This result is obviously different from the familiar Lorentz force expression 



F i = — 

1 norm . 

4"7T 



B i V j B l --V\B j B>) 



(87) 

(88) 
(89) 

(90) 
(91) 
(92) 



which follows after setting H % — B % . In normal matter we would, of course, have F^ ag = -P^orm- However, for a type II 
superconductor the combined force t p — p p V*£p due to the proton vortices must be accounted for in the magnetic force. After 
a slight rearrangement of terms we can write t p as 

1 



K = 



47ra D 



-FI 



(93) 



Comparing this result to (|9ip we see that a part of the proton vortex force exactly cancels the usual Lorentz force. That is, 
we are left with the force 



(94) 



entirely due to the proton vortices (plus a contribution in £ p due to the neutron vortices). 

In retrospect, this result is not too surprising since the global magnetic field (excluding the weak London field) is "locked' 
to the proton vortices. Compared to the standard MHD result (|92[). this effective magnetic force exhibits significant differences 
This is easy to see if we consider the special case of a non-rotating star and neglect entrainment. Then we have e ljk Vjvl = 
and it follows thai"] 

,dH cl 



9p P 



(95) 



Just like the normal matter expression (|92p . the superconducting magnetic force contains a "tension" and a "pressure" term. 
However, both terms scale with B* and H x cl — HdB 1 , instead of B z alone. Moreover, the superconducting force has an extra 
piece due to the density dependence of the energy of the proton vortices. 



This magnetic force was first derived by Easson & Pethick more than thirty years ago jEasson fc Pethicld[T977l l . for a single component 
type II superconductor using a thermodynamical approach. 
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The ratio of the two forces 



F n - 



* nor 



(96) 



will b e significant for most neutron stars. This result has been known for a long time, see for example Ea sson fc Pethick 
( 19771 ) . and several of its repercussions have been discussed in the literature. For example, as discussed bv lMen dcll (1998), the 
local speed of Alfven-type waves is modified according to v\ —¥ H c iB/Atyp p . But the overall change in the magnitude of the 
magnetic force is not the only modification imparted by superconductivity. There are also important structural differences 
between the two forces F^, s and F£ olm due to the dependence H c \ = H c i(p x ). Furthermore, we should not forget the presence 
of the force t % n — p n V l £ n in the neutron Euler equation. This force could be viewed as an effective magnetic force acting on 
the neutron fluid, as a result of the coupling due to entrainment. 

In other words, one has to bear in mind that MHD for superconducting neutron stars involves much more than a simple 
rescaling of the Alfven wave speed. 



4.5 A conservation of vorticity argument 

The variational derivation of the equations of motion is perhaps the most general and elegant treatment of the problem. After 
all, it only requires us to provide the form of the Lagrangian. Nevertheless, we wa nt to consid er the problem from a more 
intuitive (albeit less general) point of view. To do this we generalise the method of Hall (1958), originally developed in the 



context of two-fluid hydrodynamics for superfluid Helium. This approach provides a conceptually different derivation of the 
Euler equations, based on the kinematics of a conserved number of vortices carrying the entire fluid momentum vorticity. It 
also requires the input of the forces that determine the motion of a single isolated neutron or proton vortex. Thus, consistency 
between the two derivations allows us to identify the total conservative force exerted on a single neutron or proton vortex, 
without any need to study the detailed mesoscopic vortex-fluid interaction. 

The starting point of the derivation is the Onsager-Feynman conditions (|44|l and (|45[) which we can collectively write as 

W x = A/"x4 (97) 

We also need to use the fact that the vortex number density is conserved, i.e. J\f K obeys a continuity equation of the form 

d t Af x + Vj (Af x ui) = (98) 



where u x is the collective vortex velocity within a typical fluid element. In a sense, this relation defines the averaged vortex 
velocity. Taking the time derivative of (|97|) we have 

d t W x = -«£Vj (A/" x <4) + A/" x <9t4 (99) 

Reshuffling terms and using ViW x = e ljk ViV,-7r£ = we obtain 

d t wi = v, (wiui) - v, (w> x ) + af k (d tK i + <4vx - 4v 3 <) (100) 

In the language of differential geometry the motion of a single vortex is represented by the Lie-transportation of the vector 
Kx (which designates the local vortex direction) by the u x flow, i.e. 

<9t4 + £„4 = (101) 

where the Lie derivative is defined in the usual way 

£„«x = u J X7jK x — kJV jU (102) 

Then (fTOO)) reduces to 

dtWi + e y *Vj (e k i m WL<) = (103) 



which states that the canonical vorticities W x are locally conserved and advected by the u x flow. Writing this result as a total 
outer product we obtain 

a tvr j _ e ijk u*e klm V l ir? = V l * x (104) 
whe re ^ x are (so far unspecified) scal ar potentials. This result coincides with the classic formula due to Hall ( 19581 ). see 



also Siderv. Andersson fc Comer ( 20081 ) . It is important to emphasize that (I104|l applies to any condensed constituent with 



vortices, in the present case x = {n, p}. 
Alternatively, we can write (|104[) as 

o i ilk x t-tZ m r-ri ,-r, \r ilk x xv / 1 nr\ 

n x d t TT x — e J rijekimV tt x — n x V w x = A x p x e HjW k (105) 

where w xv = v x ~ u x . The functional form of the left-hand-side of this equation coincides with the vortex-free equations 
of motion (|23[) after a suitable identification of the potential \P X . The right-hand side term appears only in the presence of 
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vortices. It is not difficult t o trace its origin if we recall t hat the Magnus force exerted on a vortex (per unit length) by the 



associated fluid is given by (|Carter, Langlois fc Prix! 2001 



/l Ilk X XV / -i r\r'\ 

Mx = -Px£ KjW h (106) 

Thus, we can identify — A/x/mxi the right-hand side of (|105[) . as the averaged reaction force exerted on a fluid element by the 
vortex array. 

At this point it is clear that we can relate the macroscopic hydrodynamics we have discussed so far with the mesoscopic 
dynamics of individual vortices. This can be achieved by simply taking the force balance for a single vortex (for simplicity 
assuming a negligible vortex inertia), solving it with respect to the Magnus force, and subsequently inserting the result in 
ifTOS]) . 

We will assume that the force balance for a single neutron or proton vortex takes the form 

/Mx + /cmx + /sx = (107) 

where /c mx is a Lorentz-type force due to the interaction of the local vortex magnetic field and the charged fluids, and / s ' x 
represents "tension" force due to the self-interaction when the vortex is bent. Additional "mutual friction" forces may enter 
(|107[) . but these will be considered separately in Section [5] 

The explicit form of the tension self-force, at leading order in the vortex curvature, is 

/sx = £vxftxVjK x (108) 

where £ vx is the total energy (per unit length) of a straight vortex. 

Returning to the Euler equations (|105|) we eliminate the Magnus force using (I107[) . and expand the canonical momenta 
7r x into their hydrodynamical and electromagnetic pieces. This leads to 



n x (dt + w x V,j p x + %V ( imu x + eA - * x — p x u* j + p x ^ (£x y ^J x ) V*u x 

y 

- en x ( E* + ^ k v*B k ) = F c l mx + T x \ (x, y = {n, p}) (109) 



where we have defined the averaged vortex Lorentz force 

FLx=A/-x/Lx (110) 

and the averaged vortex tension (recall that W x = k x ) 

V=^fL -> Ti = lhLwivM (m) 

VKdp 

The explicit Euler equations for the neutron and the combined proton-electron fluids are 

n n (dt + fJV,j Pn + n„V (^rnv'i - * n - p 3 n v] ^ + p n e n w^ n \7'v J n = F c ! mll + 12 (112) 

and 

n p (d t + vpVj ) p P + n p V l (J^ mv l + 9xA) - * P - P 3 P Vj + ^ ~ /5pe P w^ n V l u P = F c l mp + T p l + ^-e 3k w^ c B k (113) 

These expressions should be compared to the equations of motion (|59[) and (|60f) (in the limit — 0), in combination with 
the eqns. (|71[) and (|80[1 for the forces t x and the potentials £ x , and the definition of the Lorentz force (|36p . This comparison 
shows that we can write 

ti = tfj - -^—e ijk tfe klm V (h vx w? - ) }■ ( 1 J n 



The term inside the bracket represents the force — /m x f° r a single vortex. Therefore this term must be equal to /c mx + /. 
as required by (|107[) . 

If, for the moment, we assume uniform densities (which also implies uniform entrainment) we have 

*x = A M T^^Wx + - c 4>Ji k tfwl e ^ = v + a; j^xe !jfc £>r \ (115) 

tain 

GTlp i ilk DC 1 ; life • "X 



and, after identifying the tension force (|108[1 . we obtain 

/cmx = ^-<l> x e vk K*wl e = --^e^jjkl (116) 
The final step we need to take is to identify 



* n = - (a*» + Cn + m$) + -mu» - p>° (117) 



*p = - (Mp + Cp + "1$) + o"™p + ej4 o - Pp u j ( 118 ) 
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The detailed identification of the various vortex forces participating in the fluid Euler equations also offers a different 
perspective on the Lorentz force cancellation result we discussed earlier. First we decompose Fl as 

Ft = Fi + ^ A/J ^ k k*j k \=Fi- £ jV«/Lx (119) 

where F£ is the Lorentz force due to the interaction between the total current and the London field alone. Then (|113p becomes 

n p (d t + VpVjJ p p + n p V (fi p + fi c + C P + m<&) - p p e p w? n V i vp = 

= f£- ^MtfL* + V+ KfL P = H + Tp - M n fL n (120) 

That is, the Lorentz force cancellation (modulo the weak London component F^) occurs due to the exact mutual cancellation 
of the terms representing the electromagnetic interaction between the proton vortex array and the charged particles. 

Eqns. (| 1 15f) and (|1 16[) provide the explicit mesoscopic forces acting on a single vortex (together with the Magnus force). 
They were "derived" under the approximation of uniform fluid densities, the same approximation that we would normally 
have made had we wished to derive these forces from first principles, that is, by a proper study of the local mesoscopic vortex 
dynamics. However, the procedure presented here allows us to move beyond the uniform density approximation and obtain the 
vortex forces to a higher precision. As before, this can be done by requiring consistency between the intuitive Hall derivation 
and the more rigorous variational approach. We can then make the following identifications for the forces acting on a single 
vortex: 

/emx = ^e'^KjtklmV 1 (0X&L ) (121) 

and 

/sx = — e" KjefcimV (£ V xK™) (122) 

At the end of the day we have derived a set of explicit expressions for the conservative forces acting on a single neutron or 
proton vortex. These results will play a key role in the derivation of the mutual friction coupling forces between the fluids. 



5 MUTUAL FRICTION COUPLING 

To complete the multifluid MHD equations, we need to account for dissipation and coupling due to the vortices. This so- 
called "mutual friction" is a key damping agent for relative fluid motion. It tends to be introduced phenomenologically 
( Alpar. Langer fc Sauls] Il984 iMendeUl Il99ll ; lAndersson. Siderv fc Comer! 120061) . This is natural since the associated micro- 



physics is complicated and not well understood. So far, we have included the mutual friction "symbolically" in the form of 
the coupling forces F^f, Fp e and FJ e (see Section f3 . 1 [I . We will now discuss phenomenological expressions for these forces and 
provide the final set of MHD equations for a mixed neutron superfluid/proton (type II) superconductor. 

The prescriptio n for deriving the mutual friction force is well established. It was first introduced in the classic work of 



Hall fc Vinenl (| 19561 ) in the context of Helium superfluidity. Initially, the mutual friction is expressed as a resistivity which 



depends on the relative motion between a given neutron or proton vortex (moving with velocities u l n and up, respectively) 
and the various fluids. By considering the equation of motion for a single vortex (essentially a statement of force balance since 
the inertia of a vortex is negligible) it is possible to express the 4 velocities in terms of the (smooth averaged) fluid velocities 
v x . These relations are then fed back into the original force expression, which now depends only on the fluid velocities « x and 
phenomenological coupling coefficients. 



5.1 The mutual friction 

We will consider two principal fric tion mechanisms. First of all, electrons can scatter dissipatively off of the vortex magnetic 



field IjAlpar. Langer fc S auls 1984). The associated drag force (per unit vortex length) takes the form 

fL = Rx(vt - <) (123) 



The coefficient R x follows from the relaxation time for electron flow through an array of vortices ( Alpar. Langer fc SaulsjB-984 ) 



and depends on the magnetic fluxes <f>*. For the neutron vortices, in particular, it will obviously depend on the strength of 
the entrainment. In both cases the resistivity (|123[) will tend to relax the vortex motion to the electron flow. 

In general, the coefficients i? x depend on the proton and neutron number densities (n p ,n n ) and will vary considerably as 
we move from the crust-core interface to the inner core. They are also sensitive (through their dependence on the size of the 
vortex cores) to the finite ambient temperature T of the neutron star core. This will be particularly important in regions near 
the superfluid phase transition, where the thermal energy is comparable to the neutron and proton pairing energy gaps, i.e. 
feflT ~ A x . These complications have been ignored in many calculations, and the drag coefficients simply taken to be uniform. 
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This simplification is not really required in the derivation of the mutual friction forces, but it makes it easier to include these 
forces in calculations of neutron star oscillation modes etcetera. 

The second mutual friction mechanism that we consider is due to the direct (magnetic) interaction between neutron 
and proton vortices. This interaction is presently only understood at a rudimentary level. Simple back of t he envelope 
estimates suggest that it coul d be strong enough to lead to neutron vortices "pinning" to the proton vortices I Sauls 19891 : 
Ruderman. Zhu fc Chenl 19981 ). If the relative neutron- proton vortex velocity is high enough to prevent pin ning, the m otion 



may still be efficiently damped as a result of the excitation of Kelvin waves propagating along the vortices ( Link 2003). 

Given the lack of a detailed neutron-proton vortex interaction force, the best we can do is adopt a phenomenological 
approach. Hence, we assume a resistive force of the usual form 

/* = C v (4-<) (124) 

where C v in another unspecified coefficient. This form for /* is undoubtedly a gross simplification. However, the subsequent 
analysis will not be straightforward if the force takes a significantly different form. If, for example, C v were to have some 
intrinsic directionality (perhaps in the form of a tensorial Cl J ) the final mutual friction forces would be significantly more 
complicated. One may, in fact, have to develop an entirely different strategy to deal with the problem. Future work should 
address this issue. We need a derivation of the neutron-proton vortex interaction based on the detailed meso-scale dynamics. 

In principle, we should also consider the back-reaction force exerted by each neutron vortex on the proton vortices. 
However, the force acting on an individual proton vortex will typically be much weaker than (|124[) . and moreover, only a 
small fraction of the proton vortex population participates in the interaction. Basically, this is due to the enormous difference 
in the density of the neutron and proton vortex arrays. The separation between individual neutron vortices is roughly 



dn ~ AC 



-1/2 



4 x 10" 



1 ms 



1/2 



where P is the spin period. Meanwhile, the distance separating proton vortices can be approximated by 



d„ ~ AC 



-1/2 



5 x 10" 



B 



10 12 G 



-1/2 



(125) 



(126) 



If we consider a fluid element of size ~Iwe would expect a number N n ~ (L/d n ) 2 and N p ~ (L/d p ) 2 of neutron and proton 



vortices, respectively, within the element. Assuming a number N p i n 
element, we see that 



L/d p of interaction sites per neutron vortex inside the 



iVpin 



2 x 10 a 



B 



1/2 



L 



10 12 G 

Comparing this to the total number of interaction sites per proton vortex 



1 cm 



N tot /N p = N pin N n /N p ~ 3 x 10" 



B 



-1/2 



10 12 G, 

we see that the backreaction onto the proton vortices can safely be ignored. 



1 ms 



L 



1 cm 



(127) 



(128) 



5.2 Vortex kinematics 



As we have already discussed in Sections 4.2-3, the motion of an individual vortex can be determined once all forces acting 
on it have been identified. In addition to the resistivity there are conservative forces present. One of these is the Magnus 
force, which is present when a vortex moves relative to a "background" flow. We also need to account for the conservative 
electromagnetic interaction between the magnetic field of a given vortex and the charged particles. Each vortex is also affected 
by a tension self-force, see Section 14.51 However, as discussed in the same Section, if we assume that the vortex arrays are 
locally straight we can extract the force acting on a single vortex from the hydrodynamical forces (|79[) simply by dividing by 
the relevant surface density Af x . This procedure leads to the vortex forces (|121[) and (|122[) . 

Having identified the forces acting on each vortex we can write down the corresponding force balance equations. In the 
case of a neutron vortex we have, 



Meanwhile, for a proton vortex we get 



fyin + /omn + fsn + /dii + /v — 



/m p + /emp + /sp + /Dp — 



(129) 



(130) 



In ord er to invert (|130[) and obtain the proton vortex velocity u p we follow the strategy from the classic work of lHall fc Vinen 
( 19561 ). Taking repeated cross products of (|130|l with rep we arrive at 



- v e + 



i + n 2 



(131) 



MHD of neutron star cores 19 



where 

fl = e iik k?wf + — (fL P + fi P ) (132) 
We have also introduced the dimensionless resistivity 

TZ P =^- (133) 

In order to obtain the neutron vortex velocity we first insert (|131[) in (| 129(1 . This leads to 

e ijh k](ul - t£) + 7L(vi -K) + D l = (134) 

where we have defined 



p n K 

and 



D l = — ( fLn + fL )+C{ 1Z p fl + e^klU } (135) 



C v = — , TZ n = — and C = - C " , TZ = TZ n +C v (136) 

p n K p n K l + Ttp 

Basically, there are three dimensionless coefficients in the problem. 
Repeating the previous inversion procedure we arrive at 

< = vl + ( nf + £lJfc ^ /fc + k^ D3fe ^j (137) 

where 

/ = e KjWk + D (138) 
According to (|131[l and (|137p we should, in general, expect relative motion between the vortices and the fluids. 

5.3 The mutual friction forces 

The final step is to translate the forces acting on a single vortex into hydrodynamical (smooth averaged) mutual friction 
forces. This is a straightforward operation if the vortices move collectively and the respective arrays are locally straight (at 
the level of each fluid element). However, this is only the simplest situation (see the footnote in Section [4T}. Nevertheless, we 
make this assumption in order to make immediate progress. It means that we simply need to multiply each force expression 
by the corresponding vortex surface density. The resulting forces (per unit volume) can be identified with the coupling forces 
in the Euler equations (|39|l . (|59|l and (|60[) . This way, we arrive at 

Fi f = Afn (/Dn + ft) = PnWn [ftn(«* - <) + C v (< - <)] (139) 
FL = Afnfhn = PnWXK - K) (140) 

F; c = AAp/i p = p P W p H p (vl - <) (141) 
In these expression we can eliminate the vortex velocities using (|131[l and (|137[l . The final results are 



= PnW n 



(142) 

/•V, = V ' Y%^( R f + e ijk tffk + ^D j kl ) (143) 



'V ^-/vWp^jl^p/I+r'"^/./ j (HI) 



6 SUMMARY: THE COMPLETE SET OF MHD EQUATIONS 

At this point we have completed our analysis of the problem of a mixed neutron superfluid and a type II superconductor. We 
have discussed the issues that arise due to the presence of vortices and magnetic fluxtubes, and accounted for the resulting 
forces in the macroscopic equations of motion. To conclude the discussion, we will now collate the information into a "complete" 
set of MHD equations. 

The final results from the previous section are rather complex. However, the model retains force terms that will be very 
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small under most circumstances. In many realistic situations the forces can be simplified considerably leading to a more 
practical MHD "toolkit". This simplified model is presented below while, for completeness, the general equations are given in 
Appendix [Bl 

In general, we have arrived at the following set of equations; First of all, provided that the constituents are individually 
conserved (i.e. ignoring reactions), three continuity equations relate the number densities n x to the transport velocities v^; 

dm* + Vi(n x <) = (145) 
for neutrons (x=n), protons (p) and electrons (e). We also have the two Eulcr equations (|59[) and (|60|) . that govern the 
evolution of the "momenta" . These Euler equations take the form 



(dt + viVi)[vl + EnW pn ] + V' (fin + $) + En^V*^ = fLi + ~ *n ~ V^n (146) 

Pa 

(8 t + «jVj)[t4 - e p wU + V 1 (p + $) - epii^vV = -^fL + fL g (147) 

Pp 

where w^ y = — v y . 

Let us now consider possible simplifying assumptions. First of all, we know that the London field vanishes identically if 
the electrons are constrained to move with the protons, i.e. when w pc — 0. In general, we see from the estimate (|58p that 
&l <C H c i , H vn and hence we can neglect the London field contribution to the forces ij. in the momentum equations in most 
cases. A similar argument shows that we can use 

W;^a p B l ->■ B l » Af p ct>oW; (148) 

where a p = e/mc and <j>o = hc/2e. Recall that we have defined 

Wl = e ijk V d K + e a viU\ = ^.W; (149) 

Wp = e ijk V 3 - [v P - e p w pn ] + a p B' = M p W l p (150) 

where A/" n and N p are the neutron and proton vortex densities (per unit area). With these approximations we obtain a 
significantly simplified set of MHD equations. First, we have 



and 

Meanwhile, the magnetic forces are 
and 



47ra D 



i B 
p 4tt 



VVSVj (H vn Wl) ~V'H V 



= — ~ — f* — VV ~ ^ 

mag — r mag ~ f,p V l^p ~ 

p P p P p 4ttp p 



47rp n 



The mutual friction follows from 



/mf 

with 

ft 



1 

Pn 



W n 



PpWp p 47ra p p p 



B J V, lH cl B' ) -\7'H cl 



I + 111 



TZpfl + e* 3fc W P /* fe 



f w e IJ ' fe W> pn +C7l p /: +Ce 1]k B ] f tk + 



47ra p p n 

Due to its relation to the current (Ampere's law), we can use the London field to infer the electron flow; 

i i 1 ijkr-j iL 

" e Vjb k 



U C ~ U P A ' 

47ra p p p 

The London field is associated with the macroscopic rotation of the two fluids, in such a way that 



6; , = -— 



1-e 



(151) 
(152) 
(153) 
(154) 

(155) 

(156) 
(157) 

(158) 
(159) 
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These relations allows us to work out the equations relating the fluid and electromagnetic degrees of freedom. They are, first 
of all, the induction equation 

1 



d t B l = e lJ e Mm Vj(v e B m ) 



(160) 



and secondly 



t^v)B h - -V 1 (fi e + m c $) + /: 



(161) 



In these equations fl represents the total mutual friction force between the electrons the other fluids. This can be approximated 
as 



p = -— If* +f* 

Pp v 



47rp p 1 + n% 



Tl P jVv^ffciB 1 ) - BV^d} - e ijk BjV k Hci + H c ie m BjB l ViB k 



(162) 



since, in a typical neutron star, J\f p 3> A/" n . 

The system of equations is closed by the Poisson equation V 2 $ = 4nGp ~ A-kG(p v + p n ) for the gravitational potential 
$, the remaining Maxwell equation ViB 1 — and by supplying the energy functional £ (the equation of state) that is needed 
to determine the chemical potentials /2 X and the entrainment parameters e x . 

These relations complete the simplified model, which should be adequate for most problems of astrophysical interest. 



7 CONCLUDING DISCUSSION 



In this paper we have developed the equations that govern the magnetohydrodynamics of superfluid and superconducting 
neutron stars. To be precise, our formalism is relevant for the conditions expected in the outer core of young and mature 
neutron stars where the matter composition is dominated by neutrons, protons and electrons. The formalism does not apply 
in the deep core, where the matter composition could be significantly different, with the likely presence of various hyperon 
species or deconfined quarks. Neither is our analysis relevant in regions of the star where proton superconductivity is not of 

type II, see the discussion in Section [2] 

In many ways, our analysis represents an upgrade of previous work on this subject, e.g. Mendel] ( 19981 ). First of all, 
we have addressed a conceptual issue regarding the dependence of the intrinsic neutron and proton vortex energies on the 
London field. Secondly, we have shown how the tension force due to the proton vortices (the fluxtubes) becomes the dominant 
magnetic force in the equations of motion while the standard Lorentz force is eliminated. We have also provided detailed 
derivations of the non-dissipative forces acting on the neutron and proton vortex arrays. Finally, we have accounted for the 
mutual friction interaction between the vortex arrays and the fluids by introducing phenomenological drag forces (in the usual 
way). 

The multifluid formalism developed in this paper is intended to provide the "next generation" of neutron star magneto- 
hydrodynamics. It can be applied to a number of interesting problems related to neutron star dynamics. Obvious possibilities 
involve the determination of equilibrium configurations of superconducting stars and the formulation of perturbation equations 
to study the oscillations of such systems. We have already discussed the properties of magnetar oscillations, and the potentially 
important effect of sup erfluidity /superconductivity on the quasi-perio dic oscillation frequencies observed during the so-called 
giant flares, elsewhere ( Andersson. Glampedakis fc Samuelsson 20091 ) . There are many other exciting problems to consider. 
Multifluid systems are known to be prone to dynamical "two- stream" instabilities provided th e re is a sufficient degree of 
relative motion between the fluid components (see, for example, Andersson. Comer fe Prix ( 20031 ); Samuelsson et al. ( 2010| )). 
Recent work for non-magnetic superfluid neutron star models has demonstrated the poten t ial importance of such instabilitie s 
for neutron star glitches and precession I Glampedakis. Andersson fc Jones 2008 fc 20091 : Glampedakis fc Andersson 20091 ) . 
Two-stream instabilities should also be present in the more realistic superfluid/superconducting neutron star model discussed 
in this paper. The nature of such instabilities need to be considered in detail. Another important application of our formalism 
would concern the secular evolution of the magnetic field in the outer core of a neutron star. Even a basic understanding of 
the magnetic field evolution in magnetar cores could tell us a lot about the way magnetic energy is distributed in the interior 
as the star ages and perhaps unveil the mechanism responsible for triggering the giant flares observed in these objects. 

Further work should also expand our formalism in two important directions. We need to consider; (i) type I supercon- 
ductivity and the transition regions between type II and type I/normal matter domains, and (ii) the magnetohydrodynamics 
of the inner core in the presence of exotic matter. These are all very challenging problems which we hope to address in the 
future. 
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APPENDIX A: INDIVIDUAL NEUTRON AND PROTON VORTICES 



In this Appendix we work out the fluid flows and magnetic fields associated with individual neutron and proton vortices, 
accounting for the entrainment in the appropriate way. We then derive the expression for the energy (per unit length) 
associated with each object. The results form a key part of the derivation of th e macros c opic f orces in Section | 4.2| 



Some of these results are well known; they can be found in, for example. iMendelll (|1991) and iRuderman. Zhu fc Chen 

1 1998^ 1. However, some aspects of our discussion have not been sufficiently exposed in the literature. In particular, we show 



that t he energy of a single vortex/fl uxtube is independent of the macroscopic London field. This result was first pointed 
out by Carter. Prix fc Langlois! ( 2000T I in the context of a simple model of a uniformly rotating superconductor with a s ingle 
fiuxtube placed along the rotation axis. A more general derivation for a mixture of superfluids was subsequently given bv lPrix 
(2000). Our analysis, which is focused on the case of a system composed of superfluid neutrons, superconducting protons and 
electrons, provides a simpler and more intuitive derivation of this result. 



Al The flow and magnetic field for an isolated neutron/proton vortex 

Let us consider the case of an single straight neutron or proton vortex. At the mesoscopic level (indicated by bars on the 
relevant quantities), we have 

^ h ViPl = <5(f), ^*Vjpg + a P K = 46(f) (Al) 

where k p = for a neutron vortex and n n = for a proton vortex. The argument of the two-dimensional delta function can 
be expressed in terms of the (cylindrical) radial distance r away from the vortex core. In our derivation we will not consider 
a detailed model for the vortex core, but focus our attention on the behaviour at distances larger than (say) the coherence 
length £ x . In essence, we consider a simple hollow-core model. 

Starting from (|A1|) we can immediately solve for the canonical momenta 

Tn = P n = o f ' n P = Pp + ~ A * = o f ( A2 ) 

where the magnetic potential is defined in the usual way, i.e. — e^VjAfc. 

In order to obtain the proton fluid momentum we first need to calculate the local magnetic field B^. To this end we need 
to consider Ampere's law 

e ijk ViBt = ^en p v p (A3) 



where the total current is supplied by the proton flow only. Taking the curl of (IA3|) and using \7jB^ — we obtain 

_j_ p,i _ 

M x A| 

where we have defined the flux parameters 



V 2 Bi-4iBi = -^6(f) (A4) 



and the London penetration length 



A^ = A« 2 (A6) 
471-pp V e / 

The entrainment parameter 6+ is defined in eqn. (|55|) . This expression leads to the estimate ([3| used in Section [2] 
The solution of (|A4)l is 

* = ^k Ko (i) ^ (A7) 

where Ko is an associated Bessel function. It is straightforward to calculate the magnetic flux through a disk surface of radius 
?~max S> A*, perpendicular with respect to k^- The relevant surface integral is 

B^dS = c/) K ki (A8) 

disk 

For a proton vortex this result represents the expected flux quantisation. 
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Going back to (|A3|) we can obtain the proton velocity (for both vortex types) 

K <fi x 



Vp 27re*A, 



In other words v p ~ exp(— r/A*) sufficiently far from the vortex axis. From the definition 

p x = m[(l — e x )^i + e x Uy] 



we find that 



Then it follows that at a distance r ^> A* 



l-£ 



Pp 



-Pn 



(A9) 

(A10) 
(All) 
(A12) 



1 - En 

for both neutron and proton vortices. 

We also need to calculate the potential A\. The symmetry of the problem dictates that = A x (r)<y3 I and we end up 
with a simple differential equation 



A. 



(A13) 



the solution of which is 



2tt 



ci 



1 



u r A 
where ci is a constant. 

For a proton vortex we obviously have p n = 0. Combining this with the previous results for Tt p and A p we find that 



(A14) 



P *= 27 



1-C! 



A 



A, 



2ttA* 



isTj ( — ) (p (proton vortex) 



(A15) 



after fixing ci = 1, in order to ensure consistency with the condition (|A12[) . Note that, we could have arrived at the same 
result without involving the vector potential by simply using eqn. (IA11|) . which for a proton vortex reduces to p p = me*v p . 



For a neutron vortex we have n p = and therefore 



Pn= 2^ 



-Al 



r - AT^ 1 U 

where (as before) we have fixed c\ = 1 in order to ensure that (|A12|) is satisfied. 



2tt l-e n 



(p 1 (neutron vortex) 



(A16) 



A2 The vortex energy 

Let us now discuss the energy £ vx of a single vortex (either associated with the neutrons or the protons) immersed in given 
macroscopic fluid flows v^. The analysis is carried out in the vortex rest-frame. In order to work out the energy associated with 
the vortex we need the difference between the localized magnetic/fluid energy and the total energies. Thus we write the total 
local momentum as the sum of p x and the (smooth) macroscopic momentum p^. We then have f4 = m x (« x — u\ + e x uiy X ), 
where u\ is the vortex velocity in the inertial frame. The total local magnetic field can be written in a similar way, as the 
sum of _B X plus the (locally smooth) London field, 

Bl oc = K+bl (A17) 
Since B x behave as ~ e~ r ^ A * sufficiently away from the vortex axis (see eqn. l|A7p ) we have the asymptotic behaviour 
B} oc -> 6L- 

The total vortex energy £ vx per unit length is the sum of its magnetic and kinetic energies plus the "condensation energy" 
associated with the normal matter in the vortex core: 

£v x — £mag + fkin + Scon (A18) 

For the present discussion it will suffice to focus on the first two terms which we write as integrals over the magnetic and 
kinetic energy densities, 

£mag = J E mag dS (A19) 

£kin = J E kin dS (A20) 

The integrals in these expressions are taken over the surface of a disk perpendicular to the vortex axis and extending to a 
radius r ~ d x , where d x is the distance to the nearest vortex. 
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For the magnetic energy we get 



That is, 



Sn 



Emas — 



B 



8 Ti 



"L 

8tt 



R 1 

8# + 4^ 



_dS + -( Ki 6 L ) 



(A21) 



(A22) 



The kinetic energy is calculated following the same prescription. Before we can do this, we need to clarify how the kinetic 
energy is d efined for a multifiuid system with constituents coupled through entrainment. Within the canonical formalism of 



Prix (2004) it can been shown that the Hamiltonian density associated with the hydrodynamical Lagrangian £h is given by 
the standard Legendre transformation Hu = n].pf — Ch- Hence, we define the kinetic energy density as 

Skm = ~(£h + Hn) = \ ( A23 ) 

x 

On the macroscopic scale, and in terms of the individual momenta, this leads to 

Skin = 2-C 1 - £ n - Ep) _1 [p n (l - £p)p n - (Pn£n + Pp£p)ppp" + Pp(l - £n)f>p ] (A24) 

In order to single out the vortex contribution we take the difference between this expression evaluated for the total momenta, 
Px + Px, an d the expression we get using only the momenta associated with the vortex itself, pi,. After some manipulations, 
where it is useful to express the London field in terms of a vector potential = e^VjA^, i.e. use 



At. = 



Pp 



it follows that 



Skin = sC 1 ~~ £n ~~ £ p) 1 [p^ 1 ~ S v)Pn + ^P( 1 _ £n ) P P _ 2 Pn£nPnP?] 



1 L» rx , Pn 



n — z 
7Pi Pn 



(A25) 



-l^ fe Vj (B^) (A26) 

47T 



When we work out the surface integral of Skin we find that the last two terms vanish due to the symmetry of the problem. 
Hence, the final expression for the kinetic energy (per unit vortex length) is, 



Pn(l — Ep) J PndS + Pp(l — En) J PpdS — 2/JnEn J PnP^dS 



4tt 



(A27) 



Each of the partial energies (|A22|) and (|A27|) show an explicit dependence on the London field. However, these London field 
terms cancel exactly in the sum £km + £mag. The final result for the total vortex energy is therefore 



£vx _ 2m2 (1 



- En 



Pn(l — Ep) J pndS + p p (l — En) J PpdS — 2p n £ n J 



plpfdS 



+ 



8tt 



BldS + £c, 



(A28) 



which is clearly independent o f the London fiel d and the macroscopic velocities. This is in agreement with the conclusions of 
Carter. Prix fc Langloisl fcoOOT l andlPrixl j2000l ). 



Inserting the results for the vortex flow and the magnetic field in the energy formula (and ignoring the condensation 
energy) we arrive at 

'A* 



^ 2 Pp 

47T£* 



lot! 



6> 



and 



47T 1 - £n 

In both cases the vortex energy is dominated by the kinetic term (the first term in (|A28|) V 
Finally, we also obtain the critical field (defined by eqn. (|78|) ~) 



47r£ vp 
n c i — — ■ 



A, 



(A29) 
(A30) 

(A31) 



Po £* 

This expression leads to the numerical estimate of Section[3] Ignoring the slow logarithmic dependence, this can be written 
as 



Scl(p P ,Pn) =h c -^-, 



h c w constant 

which is a form suitable for use in the final MHD equations of Section [6] 



(A32) 



APPENDIX B: THE FULL MHD EQUATIONS 

In this Appendix we present our final equations for the full MHD model, without the simplifications described in Section [6] 
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To begin with, and provided that constituents are individually conserved, we have the three continuity equations 

9 t n x + Vi(n x wi) = (Bl) 

for neutrons (x=n), protons (p) and electrons (e). We also have two equations that govern the evolution of the momenta. 
These Euler equations take the form 



(d t + viVj)[vt + e n w pn ] + V ! (/i n + $) + e„u4 n VX = fLf + —t' n - V*< n 

pn 

(8t + fpVj)[Up - EptllpJ + V* (A + $) - £ P WpnV'» 3 - = ~—fLf + /mag 

Pp 

where w^ y = « x — v y . 

The various force terms are given by 



47ra D 



1 - £e 



l-£n ^ 



and 



1 



/ mag = ^F mag = ^-1— J W p [v^ffdWj) - V'fflciW;)] - (V 3 6L - V l 6, L ) e* fc, V*[t;f - e pW f n ] ;> - V\„ ; L>-">) 



where a p — e/ mc and 



47ra D 



<9/9 x 9/3x 9/5x V 1 _ e " 



The canonical vorticities W x were defined as 



(B2) 
(B3) 

(B4) 



(B6) 



(B7) 



W; = e^ fc V> P - £pWpn ] + a p B l = A/pWp 

where A/" n and AT P are the neutron and proton vortex densities (per unit area). 

Recall also that the total magnetic field B % and the London field b l L are given by 



b 1 = bi + to w; - -^w n 

l 



6?. = — 



1 c n 



where O = hc/2e. 

The mutual friction follows from 



/ mf = —Fl ( = W n 



n 



with 



t P i = J^ w p [Vi(ffciWP) - V^HciWj 5 ) - (Vjtf - Vi6, L ) ] 
Using Ampere's law, we can use the London field to infer the electron flow; 

i i 1 ijk t—7 j L 

«c = " P - -. e V^&fc 

47ra p p p 

There are two more equations relating the fluid and electromagnetic degrees of freedom. The induction equation 

1 

fflp 

and 



(B8) 

(B9) 
(BIO) 

(Bll) 



(B12) 
(B13) 

(B14) 
(B15) 
(B16) 
(BIT) 
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where fl is the total mutual friction force between the electrons and the other fluids. It is given by 



Pp V 



Pn TZn 
Ppl+TZ 2 



(B18) 



The system of equations is closed by the Poisson equation V 2 <E> = AivGp w 47rG(p p + p n ) for the gravitational potential $, 
the remaining Maxwell equation Vi-B 1 — and by supplying the energy functional £ (the equation of state) that is needed to 
determine the chemical potentials /i x and the entrainment parameters e x . 
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